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ABSTRACT 

Context. In-situ spacecraft observations recently suggested that the transport of energetic parti¬ 
cles accelerated at heliospheric shocks can be anomalous, i.e. the mean square displacement can 
grow non-linearly in time. In particular, a new analysis technique has permitted the study of par¬ 
ticle transport properties from energetic particle time profiles upstream of interplanetary shocks. 
Indeed, the time/spatial power laws of the differential intensity upstream of several shocks are 
indicative of superdiffusion. 

Aims. A complete determination of the key parameters of superdiffusive transport comprises the 
power-law index, the superdiffusion coefficient, the related transition scale at which the energetic 
particle profiles turn to decay as power laws, and the energy spectral index of the shock acceler¬ 
ated particles. 

Methods. Assuming large-scale spatial homogeneity of the background plasma, the power-law 
behaviour can been derived from both a (microscopic) propagator formalism and a (macroscopic) 
fractional transport equation. We compare the two approaches and find a relation between the dif¬ 
fusion coefficients used in the two formalisms. Based on the assumption of superdiffusive trans¬ 
port, we quantitatively derive these parameters by studying energetic particle profiles observed 
by the Ulysses and Voyager 2 spacecraft upstream of shocks in the heliosphere, for which a 
superdiffusive particle transport has previously been observed. Further, we have jointly studied 
the electron energy spectra, comparing the values of the spectral indices observed with those pre¬ 
dicted by the standard diffusive shock acceleration theory and by a model based on superdiffusive 
transport. 

Results. For a number of interplanetary shocks and for the solar wind termination shock, for the 
first time we obtain the anomalous diffusion constants and the scale at which the probability of 
particle free paths changes to a power-law. The investigation of the particle energy spectra indi¬ 
cates that a shock acceleration theory based on superdiffusive transport better explains observed 
spectral index values. 
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Conclusions. This study, together with the analysis of energetic particles upstream of shock 
waves, allows us to fully determine the transport properties of accelerated particles, even in the 
case of superdiffusion. This represents a new powerful tool to understand the transport and accel¬ 
eration processes at astrophysical shocks. 

Key words. Diffusion, Shock waves, Sun: heliosphere. Acceleration of particles. Methods: data 
analysis 


1. Introduction 


The transport of energetic particles in space and astrophysical plasmas is a topic of current interest: 
for instance, the transport properties determine the measured fluxes of solar energetic particles, 


eration at astrophysical shock waves ( 

Duffv et al 

, 1995 

; Kirk et al 

1996 

Raeot and Kirk 

1997 

Perri and Zimbardo 

2012a 

Zimbardo and Perri. 

20 1?! I 

)alla et al.. 

2013 

). For many years particle 


acceleration at shocks has been described in the framework of diffusive shock acceleration (DSA). 
This theory is based on a first-order Fermi acceleration where particles, with Larmor radii larger 
than the shock thickness, can cross the shock front many times, thanks to interactions with mag- 
netic irregularities t hat scatter particles back and forth (diffusive motion) across the shock (e.g. 


Lee and Fisk. 


1982h . This process allows particles to reach high energies. However, DSA cannot 


always explain energetic particle observations in space. 

It is well known that beside normal diffusion, characterized by Gaussian statistics and by a 
mean square displacement that grows linearly in time as {x^{t)) - 2Dt, with D — \l3Av where 
A denotes the mean free path and v is the particle speed, other transport regimes are also possi¬ 
ble. In particular, they are denoted as anomalous regimes where (Ajr“) = with a I. In 

this case, the anomalous diffusion coefhcient Da has dimensions [2)„] =(Iength)^/(time)“. More 
in detail, slower processes (a < 1) are called subdiffusive, and faster processes (2 > a > 1) 
are called superdiffusive (we focus our attention on the latter). Anomalous regimes are character¬ 
ized by non-Gaussian statistics like Levy statistics, which encompasses probability distributions 


including Levy flights. Levy walks, and fractional transport equations (e.g. 

Gafter et al. 

1987 

Zaslaysky, 2002; Metzler and Klaftei, 2004; del-Castillo-Negrete et al.. 2004 

Webb et al. 

2006 

Shalchi and Weinhorst. 

2099; 

Lityinenko and Effenbergei. 

2014: 

Stern et al. 

2014 

). Numerical 


simulations of particle propagation in the presence of magnetic turbulence have indicated that 


anomalous trans 

)ort can be found, det 

tending on the particle energy and turbulence pro 

perties 

(Zimbardo, 

2005 

Zimbardo et al. 

2006 

Pommois et al. 

2007 

Shalchi and Kourakis 

2007 

Tautz 

2010 

). 


In a series of papers, 


Perri and Zimbardol (12007 


2008 


20oJbl) have shown that anomalous, 


superdiffusive transport of energetic particles accelerated at interplanetary shocks is indeed possi¬ 
ble. Those analyses were based on the observation that at few hours upstream of the shocks the 
energetic particle time profiles are power laws rather than exponentials, the latter being predicted 
by DSA. Further, the diffusion coefficient and the pitch angle scattering rate were found to be 
constant, since t he magnetic field fluctuations at the resonant scales upstream of the shocks are 


in a steady state (IPerri and Balogh , 


201C; 


Perri and Zimbardo 


2012bh . Standard arguments based 
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on DSA (IBell . 


1978h would require non-constant upstream turbulence levels to explain the ob¬ 


served power-law profiles. The latter allow the determination of the exponent of superdiffusion a; 
however, for a complete description of transport, the anomalous diffusion coefficient Da has to be 
de termined, too, as we l l as th e scale of transition to a power-law profile. The approach described 
in IZimbardo and Perril (120 13h is based on a specific type of continuous time random walk i.e. on 


Levy walks, which are random walk processes with a power-law distribution of free path lengths 
and a coupling between the free path length and travel time. This coupling is appropriate for the 
description of particles having a finite velocity. Another approach is based on fractional trans¬ 
port equations, that is on the generalization of the standard diffusion equation to the cases where 
fractional, non-integer derivatives of order jj. are present (i.e. 5,/ = K(^fld\xf^). In the usual formu¬ 
lation, fractional diffusion equations correspond to Levy flights rather than t o Levy walks, that is 
to a m odel where free path lengths and free path travel times are not coupled (IMetzler and Klaftei , 


200f)h . Therefore some differences in the prediction of the two approaches may exist. However, the 


use of fractional transport equations can be very useful, because they allow for a simple inclusion of 
external velocity and force fields, like the shock advection speed, and because powerful and well 


tested solution techniques do exist for fra ctional equations (IMetzler and Klaftei 


200C; 


Stern et al 


2014 


Litvinenko and Effenbergei 


201411 . On the other hand, the order of differentiation and the 


diffusion constant k in the fract ional diffusion equati on are often introduced “ad hoc” and thus have 
a phenomenological character (IMetzler et al.Lll998[l . Therefore, it is important to find a physically 
based way to relate atoji and Da to k. In the analysis carried out below, we show that the deter¬ 
mination of the length {q that represents the shortest path length for which a power-law probability 
distribution holds, is fundamental for characterizing the particle transport. Indeed, the values of Da 
and K are found to depend on fo- 

In this paper, we show for the first time how fo. Da, and /c can be obtained from the analysis 
of energetic particle profiles upstream of shocks. Both corotating interaction region (CIR) shocks 
observed by Ulysses and the solar wind termination shock (TS) observed by Voyager 2 are consid¬ 
ered. 

An important consequence of anomalous diffusion for particle acceleration at shocks is that 
the energy spectral index josa predicted by DSA is modified. In particular i t yield s larger spec¬ 


tral indices (softer particle spectra) in the case of subdiffusion ([K^e^^ 


1996) and smaller 


spectral indices (harder p article spectra) in case of superdiffusion (IPerri and Zimbardo 


Zimbardo and Perri. 


2012a: 


20131) . In the anomalous framework the spectral indices depend not only on 


the compression ratio of the shock r (i.e. the ratio between the plasma densities downstream and 
upstream), but also on the anomalous diffusion exponent a. We carry out a comparison between 


Perri and Zimbardo. 

2012a; 

Zimbardo and Perri. 

2013) 


2013h . We find that superdiffusion allows for a 


better explanation of the observed spectral indices. 

In Section 2 we give a description of superdiffusion in terms of Levy walks and introduce the 
main parameters that will be determined by data analysis. In Section 3 we give a description of 
anomalous transport in terms of fractional differential equations indicating the analogies between 
the two approaches and the relation between the diffusion coefficients Da and k. In Section 4 we 
present the method to extract from spacecraft data the values of a scale parameter {q, the superdif¬ 
fusion coefhcients Da and k, and the energy spectral index jssa, and we discuss the results of the 
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analysis for a number of interplanetary shocks as well as for the TS, at which energetic particle 
exhibit superdiffusion. In Section 5 we give the conclusions. 

2. Superdiffusion in terms of Levy waiks 

The interplanetary shocks studied here are large-scale shocks with a radius of curvature of about 5 
AU or more. Since the transport and the acceleration processes occur on scales much smaller than 
5 AU, these shocks are regarded as being planar, and the corresponding shock transition structure is 
considered to be one dimensional. Therefore, particles making random walks sample variations of 
the shock structure only in the direction perpendicular to the shock. To facilitate this analysis, we 
call the dimension perpendicular to the (planar) shock x and study the transport properties along 
this coordinate. 

Superdiffusion can be described microscopica l ly via a probabilistic approach that involves the 
so-called Levy random walks (e.g. 


Klafter et al 


ties of a Levy r a ndom walk, referring to 


Zimbardo et al 


(l2012h . and 


1^87 ). Here, we s ummarize the main proper - 


Metzler and Klaften (l200(lh . 


Perrone et al 


Metzler and Klaften (120041) 


120131) for more details, and to 


Zimbardo and Perri 


( 20131) for a recent analytical derivation. 


The probability i// for a rand om walker to make a free path (or jump) of length { (forward or 


backw ard) in a time t is given by (IShlesinger et al 


1987 


Klafter et al 


1993h 


1987 


Zumofen and Klaftei 


t)^A sm - vt), \(\ >4, f > 0 


( 1 ) 


where the coupling between the jump length and the travel time, as expressed by the delta function, 
is essential to ensure the finite (and, in particular, constant) particle velocity. With respect to our 
previous works, we changed the notation from // to ju -H 1. In the above expression, the power- 
law form only applies for \£\ > {q, with (q a scale parameter. An otherwise bell-shaped profile 
ensures that if/ does not diverge for £ 0. For small values of the parameter p, “heavy” tails of i/r 

are obtained, which correspond to a non-negligible probability for particles to perform very long 
jumps. It is readily verified that for // < 2 the mean square value of £, and hence the mean free 
path, is diverging i.e. {{^) - J iA((’, 0 dtdt —> oo, meaning t hat the validi^ of th e central limit 
theorem, which is pivotal in normal transport, breaks down (e.g. lKlafter et allll987h . In particular. 


for I < p < 2 we obtain superdiffusion and the parameter /r is directly related to the anomalous 
diffusion exponent a via 


a - 3 - 


( 2 ) 


Geisel et al. 

1985 

Klafter et al. 

1987; 


Zumofen and Klaftei 


1993b. 


The probability for a ‘jumping’ particle to be at a position x at time t is given by a prop¬ 
agator P{x, t). In the superdiffusive framework, this propagator is non-Gaussian, and can be 
obtained in Fourier-Laplace space by an extension of the so-called Montroll-Weiss equation 


(IMetzler and Klafterll2000il2004tlZumofen and Klaften 1 1 9931: IZimbardo and PerriLl2013h . Indeed, 


the probability distribution of the jump lengths ijjiC, t) in the Fourier-Laplace space can be written 
as 


s) = 1 - TS -I- Cl [(s -I- ikvY -H (s - ikvY], 


(3) 
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where r is the average jump time, k represents the wavenumber, and s the Laplace transform vari¬ 
able. If we make an expansion to the lowest order both in s and k, which implies the long time, 
long distance limit, and consider that 1 < // < 2, we obtain 


(^(k, s) ^ 1 - Ts H- C2\kf 


(4) 


(see 


Zimbardo and Perril (120131) for the constants Ci and C 2 ). Inserting this expression in the 


Montroll-Weiss equation and taking the inverse Laplac e transform, one obtains th e leading order 


term of the propagator for Levy walks in Fourier space (IZimbardo and Perri , 

pLw(k, t) ~ ^ exp [-Ct\kf], 

In 


201 3h . 


(5) 


where C is a constant (see below). The analytical inversion of equation (|5]) is only possible in two 
asymptotic cases. First, close to the source, that is for x one has a modified Gaussian, 




r(3/yu) 

r n 2 i 

X 

2F(l/yU) 

(cty/fl 


( 6 ) 


Zumofen and Klafterl (Il993h and 


Zimbardo and Perril (120131) for a derivation). Conversely, for 


(see 

large distances from the source, x » (Cf)'^^ but |x| < vt, one obtains to leading order a power-law. 


Plw{x, f) ^ 
(see also 


r(/t +1) 




Ct 


(7) 


Zasluvs^ 


(IBlumen et al 


(l2002h ). In addition, for |x| > vt the propagator smoothly goes to zero 


19901), this being the main difference with the Levy flight propagator given in the 


next section. 


With the additional, simplifying assumption as made bv IZimbardo and Perril (120131) . that f) 
can be approximated to be constant for \€\ < {q, it can be shown that 

1 


C = 2 




^+1 
As shown in 


cos 


(I") 


r(-//) 


( 8 ) 


Zimbardo and Perril (120131) . the superdiffusion coefhcient can be determined as 


(9) 

We can see that, basically, the value of Da depends on /r (i.e. a) and on {q. While ^ alone determines 
the exponent of superdiffusion in Eq. ([2]), both yU and the scale parameter Cq determine the values 
of the superdiffusion coefhcient that sets the scale of transport. We estimate Da directly from 
spacecraft observations. 


3. Modelling based on a fractional diffusion-advection equation 

Although the jump length and jump duration coupling introduced in the Levy walk approach above 
is in a physically strict sense more appropriate, a model based on Levy flights, which requires no 
coupling between the jump length and duration, can be regarded as an effective model for this pro¬ 
cess, similar to the diffusion approximation in classical energetic particle transport. This type of 
model has the advantage that a fractional Fokker-Planck type equation for the transport processes 
can be easily formulated, which is analogous to transport equations employed in the Gaussian 
diffusion context. This kind of equation can give solutions for the complete space-time evolution 
of the particle distribution function, which can be obtained by various complementary solution 
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techniques (ILitvinenko and EfFenbergeri l2014h . In addition, for systems where the effective ran¬ 


dom walker velocity can be different from the particle velocity, the decoupled Levy flight model 


is adequate (Idel-Castillo-Negrete et al. 


2004 


Perrone et al. 


2013 


Bovet et al. 


2014). This is the 


case, for instance, of transport perpendicular to the magnetic field due to the guiding centre motion 
caused by the E x B drift in laboratory plasmas: if the fluctuating electric field dE is strongly vary¬ 
ing and bursty, as typical of turbulence, the guiding centre velocity will be strongly varying, too, 
and the coupling between ( and t does not hold. 

Here, we only state the general idea of the model an d establish the connection to the propagator 
of the Levy walk process as d i scussed above. We refer to 


(12014 


Fichtner et al 


(12014 


M; 


Litvinenko and Effenbergen 


gdziarz^ndWeron[j2007^ 


3 E 3 


andlStern et al 


Effenbergei 

( 20lX for more 


details on the actual solutions and modelling with these equations. 

A one-dimensional example of a fractional diffusion-advection equation for shock-accelerated 
particles is given by 


d d 

—/(x, t) = t) + t) + 6{x) , 


( 10 ) 


which can be regarded as a simplified version of a more general fractional Fokker-Planck equation 


(IMetzler and Klafter 


20001). Here /(x, t) is the distribution function of energetic particles depen¬ 


dent on space x and time t, Va is a constant advection speed, and k is the constant fractional diffusion 
coefficient with the dimension (length)^/(time). In contrast to the case of normal diffusion, in the 
case of anomalous diffusion the constant k, appearing in the fractional diffusion equation, and 
appearing in the expression of the mean square deviation, are not the same. These constants even 
have different physical dimensions. Particles are injected at the shock with a delta-function source 
term. We introduced the fractional Riesz derivative with the fractional index ju as a general¬ 
ization of the one-dimensional Laplace operator. It can be expressed by the (right and left hand) 
Riemann-Liouville fractional derivatives, defined as 

1 5"' r fix') 


alfjix) = 


and 


T(m - jj) dx” 


f 


( — 1 V 


T{m - jj) dx" 




(x - x')*+^* 

fix') 


^)l+/i-, 


-dx , 


-dx ', 


( 11 ) 


( 12 ) 


with m - I < u < _ ^ where m is an int eger. The Riesz-derivative (IGorenflo et al 


Metzler and Klaftei 


2001 




1999; 


Perrone et al 


2013h is then given by 


5|x|t' 


fix) - - 


1 


2 cosifinll) 


i_^D>i + ^D>i)fix). 


(13) 


3.1. The propagator for Levy flights 

A Green’s function (or propagator) Pirix, t) for the fractional diffusion problem has to satisfy the 
fractional diffusion equation with delta source. Neglecting the advection term in Equation (fTOl) . we 
get 


^Pipix, t) = K-^^Ptrix, t) + 6ix)6it) . 
at d\xf 


Its Fourier transform is given by (e.g. 

d,PLF = -t<\krPLF + ::J-b(f), 

2n 


Metzler and Klafter 


2000I) 


(14) 


(15) 
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so that we recover the Fourier form of the //-stable distribution as 

PipiK t)^ ^ &x^{-K\kft). 

ATI 


(16) 


After computing t he inverse Fourier transform, an asymptotic expression for Plf(x, t) and x » 


is given by (ILevv 


1925 


Zaslavsky 


20021 ), 


r(// + i) 


Kt 




(17) 


thus establishing the -(1 +//) asymptotic power-law behaviour of the propagator. Comparing these 
equations with the approximate expressions for the Levy walk propagator in Equation (l5]l-(l2l), we 
see that the fractional diffusion equation can be regarded as an approximate model for Levy walks 
by identifying C as /c and having the same value of //. Therefore, the fractional diffusion coefficient 
K in Equation (fT4l i basically depends on (q and // (see Equation (O). The propagators for Levy 
flights and Levy walks differ because the tails of the latter go to zero for |x| > vt. 


3.2. The fractional superdiffusion coefficient 

Given the analogy found for the two propagators in Equations (l7]i and (fTTl i. a relation between the 
“anomalous diffusion coefficient” Da (see Equation|9]l and k can be derived. A simple dimensional 
analys is would give (Day^°' k with a - Ijp (see also the discussion ir 


Metzler and Klafter 


(I 2 OOOI) '). However, we can recover Equation |2] from the calculation of the mean square displace¬ 
ment considering a finite particle speed v. Although the formal mean square displacement diverges 
for Levy flights, one can derive a constrained or “moving” mean square displacement from the 
asymptotic expression (fTTl i for the Green’s function. Consider the definite integral 




- I x^Plf(x, f)dx 
J-vt 

2 / 7T \ 

- sin ^-//)f(1 H-//)/<•? J 


(18) 


Upon comparing with the general expression for the mean square displacement, we obtain the 
relation for the superdiffusion exponent of Levy walks, a -3— p. Eurthermore, a relation between 
the diffusion constant Da and the fractional diffusion coefficient k follows as 


Da = 


j^smg(3-a))r(4-«)|,.r 


(19) 


Using Eqs. ® and (|9]) and the reflection formula for the Euler gamma function, one can verify that 
the above coefficients are sufficiently close so that the Levy flight scenario can be used to describe 
the Levy walks to a good approximation. Thus, besides the superdiffusion coefficient, we can also 
derive k from data analysis, once having computed Da from Equation (|9l). This determination 
of K is useful for further applications of the fractional equation formalism to energetic particle 
tr ansport. The re l ation between Da and k given here is different from that shown in Equation (35) 


in 


Metzler et al. 


(Il998h or in Equation (86) in iMetzler and Klaften (120001) . which only holds for 


subdiffusion. 
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Table 1. Examples of heliospheric shocks which exhibit upstream power-law intensities of ener¬ 
getic particle fluxes 


Event s/c Date heliocentric energetic energies 

# distance [AU] particles [MeV] 

1 
2 

3 

4 


Ulysses 

11/10/1992“ 

5.2 

electrons 

0.04-0.29 

Ulysses 

22/01/1993“ 

5.0 

electrons 

0.06-0.18 

Ulysses 

10/05/1993“ 

4.7 

electrons 

0.04-0.29 

VnyagPT ? 


7 

inns, , 

-0.54-3.50 


a :IPerri and Zimbardol ll2009atl : b :IPerri and Zimbardol 120081) : 
c : 


Perri and Zimbardol (120071 ): d : lPerri and Zimbardol d2009bh . 


Table 2. Measured parameters of the heliospheric shocks listed in Table [T] 


# 

[km/s] 

[km/s] 

n\ 

[cm“^] 

112 

[cm-^] 

ObS 

s/c 

^break 

[hours] 

1 

755 ± 10 

643 ±4 

0.13681 ±0.00007 

0.3617 ± 0.0006 

68 ° ± 11 ° 

8 

2 

760 ± 1 

583 ±4 

0.110236 ±0.000004 

0.572 ± 0.006 

35° ± 13° 

20 

3 

814 ±50 

645 ±311 

0.1171 ±0.0003 

0.352 ± 0.002 

59° ± 10° 

4 

46 

300 

155 

0.001 

0.00158 

74°-82° 

240 

4“ 

368 

145 ± 112 

9.6194 ± 0.0006(10--*) 

2.0757 ± 0.0002(10-3) 

74°-82° 

240 


For events 1-3: data from CDAWeb; for event 4: data from space.mit.edu/pub/plasma/vgr/v2/dailv; 


a : from 


Balogh et all (Il99.5h . b: plasma values from 


Richardson et al 


20031 : 


c: plasma values computed over a sliding window of 7 days over a 100 days period. 


For event 4, when the statistical error is > 100% only the error values are reported. 


4. Application to spacecraft measurements 


We now directly estimate the superdiffusion coefficient 2)„ and the fractional diffusion coefficient 
K, derived from the two approaches described in Sections 2 and 3, from spacecraft observations for 
a number of interplanetary shock waves. We selected those ev en ts for which part i cle su perdiffusion 

. LoOOafl 


has already been observed (iPerri and Zimbardo 


2008 


Zimbardo et al 


\2Qm . 


4.1. The scale parameter 

Recall that the paramater fo indicates the scale at which the jump probability f) turns out to 
have a power-law shape. In the observations, this is related to the distance from the shock at which 
the particle time profiles begin to exhibit a power-law behaviour, too. There is indeed a break 
between the shape of the particle time profile in the vicinity of the shock and its shape at few hours 
upstream from the shock front (see Figure 1) and this corresponds to the two asymptotic cases for 
the propagator which have been shown in Eqs. ® and (|7]i. 

Considering the form of C in Eg. (Ell we can defin e a scaling variable ^ = (Jc/^o)/(f/fo)*^^ 


dZumofen and Klafter 


1993 


Zimbardo and Perri. 


20131), where fo = ^o/v, so that we have ^ 1 


close to the shock. The modified Gaussian propagator in Eg. ® can be used, giving a nearly flat 
profile (indeed, for ^ 1 the Gaussian is slowly varying). Conversely, f » 1 represents a region 

farther away from the shock and the power-law f orm of the propagator in Eg. dill can be used, 
giving a power-law particle profile as described bv IPerri and Zimbardol (120071 l2008h . Clearly, the 
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1993 May 10, 19:17:00 UT 



Fig. 1. Hourly resolution electron flux measured by the Ulysses spacecraft upstream of the CIR 
shock crossing of May 10, 1993 (solid line). The time Af represents the difference between the 
upstream measurement time and the shock crossing time. The power-law fit is represented by the 
thin solid line, and the time of the break in the power-law by the dashed red line. 


change from one form to the other corresponds to the case ^ » 1. Therefore, it is possible to ob¬ 
tain the distance Xbreak of the break, which separates the different forms of the upstream energetic 
particle profiles from the condition ^ = 1. Thus, we find 

.^break _/ ^break 

7o \ U 

Further, the relation Iq - vto is implied from the space-time coupling present in Eq. O, so that 
giving the scale length 7o is equivalent, for particles of given velocity known from spacecraft mea¬ 
surements, to giving the scale time fo- Since we are deriving the parameters for energetic particles 
accelerated at interplanetary shocks, we denote the upstream and downstream values by the sub¬ 
scripts 1 and 2, respectively, and the reference frame i.e. spacecraft or shock, by the superscripts 
‘s/c’ or ‘sh’, respectively. Assuming a steady-state shock, the particle motion upstream of the shock 
can be described as a competition between the supersonic plasma advective motion, x = Vf't, and 
the superdiffusive random motion (Ax^) = 22 )q,U. Then, we insert in Eq. (l20l i an expression for 
.^bieak that comes from the advective motion, namely 

.^break ~ Ej tbreak ■ (71) 



Using Eq. (l2Tb and 7o = vfo and rewriting in terms of distances, we obtain from Eq. (l20l i 

Vl/iZ-ff) 

1 / 

.^break “ ^0 


T7Sh 

V kj / 


( 22 ) 


which gives the distance of the break from the shock (see Eigure[T]). Equation (l22]i compares with 
the expression for the exponentiation length Lprec = Dl'Vf' = which is obtained in the 

case of normal diffusion for a constant diffusion coefficient (e.g. 


Lee and Eiskl 


19821). However, 


it is worth remembering that 7o does not have the same meaning as the mean free path A, which 
diverges in the superdiffusive case. The higher the upstream plasma flow velocity Vj*', the shorter 
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either Lprec or ;cbreak- This is because the advective motion is ‘squeezing’ energetic particles towards 
the shock. We may notice that studying the precursor, i.e. the prohle of upstream particles, is the 
best and most direct way to study the transport properties perpendicular to the shock (provided the 
upstream is known). 


4.2. The energy spectral index of particles accelerated at the shocks 

The transport of energetic particles also influences the shape of the energy spectrum. In partic¬ 
ular, DSA predicts an en ergetic particle spectrum that only depends on the compression ratio of 


the shock (iDrurv 


19831). This holds for both relativistic and non-relativistic particles. Thus, the 


investigation of the energy spectra for particles accelerated at shock waves can be a test for the 
validity of DSA. Recent observations during the TS crossings by the two Voyager spacecraft gave 
indications of a compression ratio at the shock of r ~ 2, although large-scale tempo ral variations 
of the TS makes it difficult to accurately estimate the value of the compression ratio. 


Decker et al 


(120081) have shown that the slope of the differential spectrum over a broad range of energy channels 
is ~ 1.26. This value can be recovered by DSA, assuming r ~ 3, which is in strong disag reement 
with the value obtained from Voyager observations. Recently, 


Arthur and Le RouxI (120131) numeri 


cally integrated the focussed transport equation for a spherical, stationary shock that includes both 
a precursor and a subshock. They assumed a global compression ratio of about 2.8 . The energy 
spectra obtained were in good agreement with those derived from the Voyager data, thus indicating 
that particles with diffusion length larger than the whole shock size (including the precursor and 
the subshock) will experience a compression ratio greater than 2. This result is consistent with a 
DSA scenario, how ever the effective compression ratio sampled by particles can be energy depen¬ 


dent ( Amato 


2014h so that low energy particles will experience a compression ratio closer to that 
observed at the subshock i.e. r ~2. 

We propose an alternative explanation based on superdiffusion for the spectral index values 
determined for particles accelerated at the TS, and for differential e nergy spectra of particles ac 


Perri and Zimbardol (l2012ah 


celerated at a couple of CIRs reported in Table [T] (see details below), 
derived the particle spectral index in the framework of the superdiffusive shock acceleration (SSA); 
in this case the spectral index depends on both the compression ratio and the index of superdiffusion 
a. For ultra-relativistic particles the differential energy spectral index is 
6 2 - a 


r = 


+ 1 , 


r - 1 3 - a 
while for non-relativistic particles 
3 2-a 


Perri and Zimbardol (l2012al) obtained 


7 '■ 


1 3 


-H 1. 


(23) 


(24) 


Perri and Zimbardol (12012ar) further calculated the spectral index yj - y-O.S of the differential flux 
djIdE for non-relativistic particles ass uming superdiffusive transpo rt; with a compression ratio 
r - 2 and a - 1.3, the value found in IPerri and Zimbardol (l2009bl) for the particles accelerated 


upstream of the TS-see Table [3] we obtained a spectral index of ~ 1.7 (see Table |4]i. This value 
is very close to the 1.67 value foun d in the low energy charged particle spectrum detected by 
Voyager 1 just after the TS crossing ((Decker et al.L l2005l). but not in a very g ood agreement with 


the 1.26 value reported from the Vo yaser 2 observations ((Decker et al 


compression ratio, as considered by lArthur and Le RouxI (120 13h . can bring the predictions closer 


20081) . Thus, an increased 
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to the observations. Still, SSA is able to give harder spectral indices so that r - 2.4 would be 
sufficient to explain the observed yj, assuming a = 1.3. Below we will compare the spectral indices 
predicted by SSA with the indices directly obtained from the observed energy spectra of particles 
accelerated at ClRs. Electrons accelerated at CIR events in Table [1] are relativistic, so that the 
measured differential flux (assuming that the particle speed within each energy channel is almost 
constant) can be written as 


di dN dN dp dp 

— = V— s V - oc vp —, 

dE dE dp dE dE 


(25) 


where jp is the index of the differential spectrum in momentum. After manipulating Eq. (l25l) . being 
E = yj it can be easily derived 


djIdE ~ p-^” 

jp = [(6/(r - 1))(2 - a)/(3 - a)] + 1. (26) 

Since for the events in Table[T]both r and a are known, yp can be computed and compared with the 
slope of the djjdE spectra observed. 


4.3. Data analysis and method 


Eor interplanetary shocks measurements are made by spacecraft as a function of time, so that we 
need to estimate Vbreak by the observed break time as seen by the spacecraft. This is not the 
break time for particle advection introduced in Eq. (|2T]i, but the break in the particle spatial profile 
as seen by the spacecraft when sampling different distances from the shock, as shown in Eigure[T] 
Thus, the break position is obtained as 


•'break - 


(27) 


where is the relative velocity between the spacecraft and the shock. 

To estimate the upstream Vj** in the shock frame we have to take the following transformation 
of velocities into account: 


ys/c ^ ya/c ys, ^ 


/Sh 


(28) 


We make the simplifying assumption that the shocks are planar and perpendicular to the radial 
direction, even though it is well kno wn that CIR shocks are oblique with respect to the radial direc¬ 
tion (e.g. lGosling and Pizzolll999h . A study of the full set of Rankine-Hugoniot jump conditions 
is deferred to future work. Considering further that the compressi on ratio is to be expressed in the 


shock frame as r = we can easily obtain the relation (e.g. [Burgess, 


199.5h 


yS/c 


rV: 


■s/c 


v; 


s/c 


r- 1 


(29) 


The compression ratio can be obtained from the observed densities as r = n 2 /ni, although strong 
fluctuations of th e measured p a ramet ers, like plasma density, v elocity, etc. are usually found (see 
the discussion in iBalogh et al.l (Il995h and in iGiacalonel (l2012h ). Eor example, in the hourly aver¬ 


aged Ulysses data, the plasma density can change by 50% in a few hours. Therefore, to estimate of 
the plasma parameters upstream and downstream of the shock waves, we calculated both plasma 
speed and density within a time window of seven hours for those CIRs detected by Ulysses for 
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which superdiffusion of energetic electrons was observed. These Ulysses events are listed in Table 
[1] After a first iteration, the window is moved by one hr, then the calculation of velocity and density 
is repeated. This process continues by spanning a 12-24 hours time interval upstream and down¬ 
stream (that is shifting the time window for 12-24 hours), depending on the “regularity” of the 
velocity and density profiles in each event. To avoid sharp variations of the plasma parameters 
close to the shock fronts, we started the above evaluation procedure after a couple of hours from 
the time of the shock, both upstream and downstream. Finally, we calculated an average value for 
the plasma speed and the density from the values computed within the running windows, along 
with a standard deviation (statistical error); these are reported in Table |2] Although the procedure 
described here for the plasma parameters calculation is accurate, the intrinsic strong variability of 
the time series upstream and downstream of the shocks can lead to large unavoidable errors in their 
estimation. When the standard deviations are larger than the averages, only the error values are 
given. 

For the TS crossing (event 4 in Table [1]), the analysis required more attention. For a fir st esti¬ 
mation, we considered the shock and plasma parameters given by iRichardson et al.l (l2008h at the 


so-called TS-3 event (i.e. shock speed in the spacecraft frame -68 km/s, and compression 

ratio of the shock r ~ 1 .58). In particular, we e stimate velocities and de nsities upstream and down¬ 
stream of the TS from the profiles in Figure 3 in iRichardson et al.l (120081) in the vicinity of the shock 


crossing. For a second estimation, we also considered a larger scale tr end of the plasma parameter s 
to be consistent with the analysis of energetic proton time profiles in IPerri and Zimbardol (l2009bl) . 
Indeed, a power-law decay for the proton time profiles was found for about 100 days upstream of 
the TS. Thus, we have further calculated the plasma densities and speeds over sliding windows of 
seven days for a 100 days period upstream and downstream of the TS crossing on 2007 September 
1st. In this case, we assigned to the shock a speed (in the spacecraft frame) of -27 km/s. This 
double estimation is because, while for TS-3 Voyager 2 crossed the inwards moving TS, which 


was a typical quasi-perpendicular shock with shock normal angle Osn = 74.3° (IRichardson et al 


20081), other TS crossings exhibited different parameters. Therefore, to estimate an average shock 
sp eed, we also considere d the results of a global MHD simulation of the TS dynamics, performed 
by[ 


Washimi et al. 


(120071) : considering the solar wind parameters of 2006 and 2007, they found that 
the shock is moving towards the sun (inwards). Thus, the relative velocity between the spacecraft 

s/c 

and the shock is given in the reference frame of the spacecraft by - IVshl + Vs/c ~ 27 km/s. 


being 14h-12 km/s the TS’s inwards speed with respect to the Sun (IWashimi et al.l l2007h . and 

Vs/c ~ 15 km/s the spacecraft’s speed with respect to the Sun. From Table [3] it can be seen that 
both the break distance Xbreak and the scale length {q depend on the assumed shock speed. 

For each particle species, i.e. electrons for events 1-3 and protons for event 4, we can obtain 
the speed from the energy 


1 


(30) 


£2 

where E = mc^ + Eh is the total free particle energy, with E^ the ‘kinetic’ energy determined by the 
spacecraft instruments. All energies can be conveniently expressed in keV. 

An expression for the scale length £o can now be obtained by inverting the relation (l22l) above, 

' ysh \ ^ / (2~“) 


^0 — -^break 


(31) 
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The value of a - 3 — u is obtained by the slope of the power-law upstream time profile for each 
energy channel ( Perri and Zimbardo . 2007 , 2008 . 2009 Jb ). To estimate the particle speed, we used 


the average energy of each energy channel. Once the value of fo is known, the values of and /c 
can be obtained. 


We further computed the particle energy spectral indices for the CIR events 1 (both upstream 
and downstream) and 2 (upstream). The downstream fluxes of event 2 and the event 3 in Table 
[1] have not been analyzed for obtaining jp (see eqs|^ l26l l since the electron time profiles fall 
down rapidly, thus making unreliable the computation of an average flux over a time window of 
few hours. Moreover, the data within the different energy channels tend to be almost overlapped, 
making it impossible to obtain a clear scaling in the energy spectrum. The reasons for the high 
variability and the overlapping of the particle time profiles are not clear and require further inves¬ 
tigation. For electron time profiles of event 1 (Figure |2]) and for upstream time profiles of event 
2 (Figure [3]), we calculated an average differential flux < djldE > over the time windows delim¬ 
ited by the vertical red dashed lines in Figures |2] and |3] (the shock positions are indicated by the 
thick vertical arrows). Figure |4] shows the obtained < dj/dE > as a function of the “momentum” 
cp = yjE^ + Imc^Ek for event 1. The kinetic energy Ek corresponds to the ‘centre’ value of each 
energy channel. The same plot is shown in Figure|5j but for event 2. The slopes of the differential 
energy spectra are reported in Table|4]along with the theoretical values expected by DSA and SSA 
once the compression ratio of the shock and the exponent a of superdiffusion are known. 


4.4. Results of data analysis 


The results of our analysis are summarized in Table |2] (measured parameters with errors) and Table 
[^(derived parameters with errors) for both the interplanetary shock crossings by Ulysses (events 1- 
3) and for the TS crossing by Voyager 2 (event 4). 


Table|3]shows that for event 1 the obtained values of fo range from 30 m to 600 km, depending 
on the energy channel. Although the particle velocity enters the expression of {q, most of this vari¬ 
ability is due to the variability of a, which determines the exponent in Eq. OTT) . and to the small 
value of the ratio Vf^jv ~ 10“^. The values of €o reported are actually those associated with sta¬ 
tistical errors, since errors are always larger than 100%. The analysis shown in IPerri and Zimbardo 


(l2009ah clearly indicates that the energetic particle time profiles upstream of the shock decay as 
power laws, rather than as simple exponential (as envisaged by DSA), thus superdiffusion is ex¬ 
pected. However, the exponent of the mean square displacement a, which gives the degree of 
superdiffusion, can vary according to the particle energy. Indeed, different energies imply different 
Larmor radii and then a different (resonant) interaction between particles and magnetic fluctua¬ 
tions. In turn, this implies the large variability of {q. The values of the diffusion coefficients Da and 
K (computed using the (q values in Table [3]) change in the same “direction” as Iq, but with much 
smaller deviation. In particular, the values of Da are rather stable, a property that is important for 
developing predictive algorithms for the fluxes of solar energetic particles. 


A similar behaviour is found for events 2 and 3. On the other hand, event 2 is only slightly 
superdiffusive (a = 1.0-1.19) and this gives larger values of io and of the diffusion coefficients. In 
the limit of normal diffusion. Da and /c cannot be computed because r(a - 1) diverges. While this 
divergence hints at making the above coefficients closer to the large diffusion coefficient usually 
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Table 3. Derived parameters for the heliospheric shocks listed in Table [T] 


# 

yS/c 

*^sh 

ysh 

[km/s] 

Vf 

^break 

[10®km] 

Ek 

[keV] 

V 

[“] 

a 

4 

[km] 

Da 

["] 

K 

n 

1 

576 

179 

68 

16.6 

53.5 

127 

1.44 ±0.08 

620 

450 

4 






88.5 

157 

1.56 ±0.09 

30 

220 

0.5 






145 

188 

1.7 ±0.1 

0.03 

80 

0.04 






234 

218 

1.6 ±0.2 

7 

300 

0.3 

2 

541 

219 

42 

38.9 

53.5 

127 

1.00 

_ 

_ 

_ 






88.5 

157 

1.08 ±0.02 

15000 

20000 

10000 






145 

188 

1.19 ±0.03 

9000 

5000 

600 






234 

218 

1.02 ±0.05 

30000 

130000 

100000 

3 

560 

219 

73 

8 

53.5 

127 

1.29 ±0.08 

9000 

2000 

100 






88.5 

157 

1.38 ±0.07 

2000 

1000 

20 






145 

188 

1.31 ±0.08 

4000 

2000 

100 






234 

218 

1.15 ±0.08 

15000 

10000 

2000 

4 

-68 

368 

233 

59 

765 

12 

1.30 ±0.07 

400000 

1500 

20 






1565 

17 

1.29 ±0.08 

270000 

2000 

25 






2820 

23 

1.3 ±0.2 

100000 

1500 

15 


-27 

396 

183 

23 

765 

12 

1.30 ±0.07 

180000 

3500 

300 






1565 

17 

1.29 ±0.08 

100000 

4500 

350 






2820 

23 

1.3 ±0.2 

60000 

2500 

150 

Notes on units: a : 

10^ km/s; b : 10*’ km^/s“; c : 10® km''/s. Errors have also been estimated for derived 


quantities (not shown). 

For fo, Da, and k the statistical errors are larger than their mean values, therefore only the error values are 

reported. 



days of 1992 


Fig. 2. Hourly resolution electron flux measured by the Ulysses spacecraft at the CIR shock of 
11/10/1992. The shock position is indicated by the vertical thick arrow. The vertical red dashed 
lines delimit the regions upstream and downstream where an average flux < djldE > has been 
computed. 
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Fig. 3. Same as Figure|2] except for the CIR shock of 22/01/1993. 


obtained for normal diffusion, one has to remember that most of the expressions given in Section 
2 are valid for a strictly larger than one. Also, event 2 is only marginally superdiffusive. This can 
be explained by considering that this event is a quasi-parallel shock (see Table |2]- d^n ~ 35°). 
Thus, an upstream ion foreshock is formed, which enhances the level of turbulent fluctuations and 
wave activity due to the backstreaming reflected ions. This results in a stronger interaction between 
particles and fluctuations so th at there is a higher probability of recovering a normal diffusive 


scenario (IPommois et al. 


20071). 


In the case of the TS event, a shows similar values in all the energy channels; in this case fo, 
Da, and k do not broadly vary but tend to be more stable. 

Despite the large variability of {q, we note that Xbreak, which is directly determined from the 
spacecraft data, does not exhibit such a large variability but is quite stable from shock to shock. 
It is important to remark that {q is the particle jump length above which the probability distribu¬ 
tion of the particle jump lengths i.e. f), decays as a power-law (see Eq.(3)). For each event, 
Table[3]shows that higher values of the exponent a of the mean square displacement correspond to 
smaller values of {q. When the value of a increases (but always bounded 1 < a < 2 in superdiffu¬ 
sive processes) the particle propagation becomes more and more superdiffusive, therefore £o might 
become shorter and shorter, so that the particle during its propagation has a higher probability of 
making both very long and very short jumps. This gives rise to a very non-homogeneous, scale-free 
propagation of particles in space. 

The electron energy spectra for the CIR of 11/10/1992 are shown in Figured in particular, the 
top panel in Figure|4]displays the energy spectrum (filled circles) derived from < djIdE > upstream 
of the shock, while the bottom panel reports the spectrum computed downstream. The best power- 
law fits are indicated by the thick solid lines (the slopes are also displayed in the panels). The 
black dashed lines indicate the trend expected from DSA, while the red dotted-dashed lines refer 
to the power-law behaviour predicted by SSA considering the minimum and the maximum a value 
found by the analysis of the electron time profiles (see Table |3]l. It can be noted that the observed 
energy spectra, both upstream and downstream, are well described by the SSA model, although 
the agreement is much better upstream. This can be explained by the fact that we have calculated 
jp downstream using the values of a found for particle time profiles upstream. The same plot has 
been produced for the CIR event of 22/01/1993 (see Figure [S]). However, in this case neither the 
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11 /10/1992- downstream 



Fig. 4. Energy spectra (filled circles) for the CIR shock of 11/10/1992 calculated upstream (top 
panel) and downstream (bottom panel). The thick solid line represents the best power-law fit. The 
black dashed lines indicate the DS A prediction and the red dash-dotted lines indicate the power-law 
behaviour predicted by SSA considering the minimum and the maximum values of a found from 
the electron time profiles (see text). 


DSA and the SSA are able to reproduce the observed electron energy spectrum, which results to 
be much steeper. Notice that for this event the compression ratio is > 4 (see Table |4]i, suggesting 
non-stationarity and a large time variability for the plasma parameters upstream and downstream 
of the shock . As d iscussed above, a comparison between the Voyager observations reported by 
Decker et al.l (l2008h and the values predicted by the theoretical models gives indication of a better 


agreement with the SSA model (as indicated in TableHJi. 
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Fig. 5. Same as Figure|4] except for the CIR shock of 22/01/1993. 
Table 4. Results from the analysis of the electron spectra. 


# 

r 

a 

Tobs 

7dsa 

Tssa 

1 u 

2.6 

[1.44-1.70] 

2.26 

2.9 

[1.87-2.35] 

1 D 

2.6 

[1.44-1.70] 

1.43 

2.9 

[1.87-2.35] 

2U 

5.19 

[1.02-1.19] 

3.7 

1.72 

[1.64-1.71] 

4 

2.0 

1.3 

1.26 

2.0 

1.7 


5. Conclusions 


In conclusion, we have presented the description of anomalous, superdiffusive transport in terms of 
a probabilistic description, corresponding to Levy walks, and in terms of a fractional transport equa¬ 
tion, corresponding to Levy flights. Although they represent two different approaches to anomalous 
transport, we have shown similarities between them and also how the mean square deviation fore¬ 
seen for Levy walks can be recovered from the solution of the fractional transport equation. These 
results are interesting because the use of the fractional transport equ ations has extensively been 


applied to a variety of physical systems, including laborato ry plasmas (Idel-Castillo-Negrete et al 


2004 


Bovet et al. 


2014 


Litvinenko and Effenbergen.l2014t) . Further, it is important to understand 


to what extent Levy flights can be used to recover the properties of Levy walks. Within those frame¬ 
works, parameters of superdiffusive transport, as the corresponding diffusion coefficients Da and 
K, have been expressed in terms of the scale parameter {q. The latter and, in turn, the superdiffusion 
and the fractional diffusion coefficients have, for the first time, been estimated directly from space¬ 
craft observations at heliospheric shocks. The values of Iq reported in Table[3]are fou nd to be much 
shorte r than the mean free path A which is, according to the standard estimates (e.g. 


Bieber et al 


19941) . of the or der of 0.1-1 AU a t a heliographic distance of 1 AU, and of the order of a few tens 
of AU at the TS (lGiacalonell2013h . These low values of the scale length fn have been in terpreted as 


Zimbardo and Perril (l2013h l. It is how 


indicative of very short acceleration times (see Eq. (40) in 
ever important to recall that {q does not have a meaning similar to the mean free path A: firstly, the 
mean free path is diverging in the case of superdiffusion, and secondly, {q is the shortest length for 
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free paths having a power-law distribution. This means that free paths much longer than {q have 
relative high probability, and their occurrence gives rise to superdiffusion. We have also checked 
the expression of the spectral index predicted by SSA for the considered shocks. Within the exper¬ 
imental uncertainties, SSA gives a much better explanation for the observed yobs values than DSA, 
thanks to the capability to give harder spectral indices than DSA. The results obtained give insight 
into the random walk properties for superdiffusion, furnish an estimation of the input parameters 
for numerical simulations, and supply the basis for further studies of energetic particle acceleration 
at heliospheric and astrophysical shocks. 
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